clear
clc

snr='none';
ibo=10; %The input power back off
loop=10000; % iteration times

me=modem.qammod(16); %modulation scheme
de=modem.qamdemod(16); %demodulation scheme
N=2048;  % OFDM subcarriers number
cap=N;   % Look up table size
num=sqrt((10)*10^(ibo/10)); %calculate the power back off
aa=0.1; % iteration steps

Tab_a=ones(1,cap); % look up table for amplitude
Tab_p=zeros(1,cap); % look up table for phase
address=zeros(1,N); %  address table
e=zeros(1,loop); % error for each iteration
%--------------------------------------------for speed :)
xTa=ones(1,N);
xTp=zeros(1,N);
Ina=ones(1,N);
Ina1=ones(1,N);
Inp=zeros(1,N);
Outa=ones(1,N);
Outp=zeros(1,N);
%-----------------------------------------------
for start=1:loop
  
    data=randint(1,N,16);
    data1=modulate(me,data); %generate 16QAM signals
    xo=data1/num;  %signal power back off
    d=ifft(xo)*sqrt(N); % OFDM signals
    Ina=abs(d); 
    Inp=angle(d);
    for uj=1:length(Ina)                          % time zone predistorter
        [address(uj),Ina1(uj)]=Q(Ina(uj),1,cap);% Q: adrress function
        Outa(uj)=Ina1(uj)*Tab_a(address(uj));
        Outp(uj)=Inp(uj)+Tab_p(address(uj));
    end

    [ya,yp]=HPA(Outa,Outp);

    y=ya.*exp(1i*yp);

    for ui=1:length(address)    % the look up table renew algorithm
        Tab_a(address(ui))=Tab_a(address(ui))-aa*(ya(ui)-Ina1(ui));
        Tab_p(address(ui))=Tab_p(address(ui))-aa*(yp(ui)-Inp(ui));
    end

    %----------------------------------------------- end of transmitter

    %------------------------------------------------------------channel
    if strcmp(snr,'none')
        Yreceive=y;
    else
        Yreceive=awgn(y,snr,'measured'); % add noise
    end
    %-------------------------------------------------------------end channel
    %------------------------------------------------receiver
    Go=Yreceive;     % feedback

    Gone=fft(Go)/sqrt(N);  %  demodulation of OFDM
    Gone=reshape(Gone,1,length(Gone));  

    Result=demodulate(de,Gone*num); % demodulate
    z=Gone*num; % the signal constellation
    e(start)=mean(abs(Gone-xo).^2);
    EVM(start)=mean(abs(d-y))./mean(abs(d));
    [bes,ber(start)]=biterr(Result,data);
end
BER=mean(ber(loop-2500:loop));
figure(1)
plot(z,'*') %constellation figure
